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ABSTRACT 

For nearly two decades we have witnessed an intensive development 
of a statistical methodology for assessing length of life and relia- 
bility of performance from empirical data. The initial stimulus for 
research on statistical problems in life testing and reliability came 
from the need to answer pressing practical questions which could not be 
treated by the existing statistical techniques. Because life and per- 
forraance tests are so time consuroing and expensive to run^ it is a 
practical necessity to terminate them as soon as possible. 

For the statistician this means developing estimation and decision 
procedure for data, which are severely curtailed in one way or another 
long before all items on test have actually failed. The 
estimation is more complicated when the data are truncated, i.e. when 
the observer loses track of some individuals before death occur. The 
product limit method of Kaplan and Meier is one way of estimating p(t) 
when the mechanism causing truncation is independent of the mechanism 
causing death. 

This paper proposes alternative estimators and compares them to 
the product limit method. A computer simulation is used to generate 
the times of death and truncation from a variety of assumed distribu- 
tions. No single estimator gives the best fit to the "true'* distribu- 
tion of death under all situations. However, other estimators are 
shown to be better than the product limit estimator under all of the 
assumed situations. 
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I. INTRODUCTION 



Let the random variable T denote the time that elapses until an 
event occurs; the event may for example be an equipment failure, an 
individual's death, or the detection of a target. Denote by p(t) the 
probability of survival to time t, 

p{T > t} = p(t) 

Picturesquely, T is called a lifetime, and p(t) is a survival pro- 
bability; 

F(t) = 1 - p(t) is the distribution function of T. 

In the medical field, one might wish to estimate the probability, 
p(t) that a patient survives t after a certain surgical procedure 
for cancer. In electronics, one wishes to estimate the probability of 
continuous failure-free operation of an equipment for time t. In the 
military, one might be interested in the probability of conducting a 
certain mission, under specified environmental conditions, without de- 
tection by the enemy. The event of interest may be a human death, 
equipment malfunctions, or sonar detection. Following Kaplan and Meier, 
Reference (1) , this paper will refer to the event of interest as a 
"death". The test element in the sample may be a human, a radio, or 
a submarine. This paper will refer to the test elements as "individuals". 
Suppose that observed values of T are tw t_, t_,...t , so that N 
lifetimes are observed. In this case an appropriate (unbiased) esti- 
mates of survival to time t is 



6 



number of t. *s > t 
p(t) i 

N 

Under many circumstances complete lifetimes are not observed; censoring 
occurs at certains, , beyond which the life of an individual is not 
known. In such cases construction of an appropriate estimate of the 
survival probability is more difficult. In this paper various estimates 
of survival probability are studied when lifetimes are randomly censored. 
This means that censoring times are assumed to be realizations random 
variables independent of the actual lifetimes. 

The product-limit estimator of Kaplan and Meier, Reference (1), is 
an accepted method of dealing with the problem of censored data. This 
paper presents thirteen non-par ame trie estimators, including the product 
limit function. Censored data sets are simulated. The thirteen esti- 
mators are compared by examining their performance on the simulated data 
bases . 
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II. THEORY 



There are two approaches to the empirical estimation of the survival 
probability/ p(t): 

(1) one may use the observed fraction of survivors at arbitrarily 
selected times (step function estimator) / or 

(2) one may focus attention on the times of the observed deaths 
(point estimator). 

The initial discussion is based on the assumption that all observa- 
tions are complete, i.e., it is assumed that all individuals remain under 
observation until their time of death. This initial assumption is for 
the purpose of simplifying the discussion. Then, later in this paper, 
the discussion is broadened to include incomplete data with observations 
of both death and censoring events. 

Survival Probabilities; No Censoring 

Let 0=t^ <t, <t^ ... <t. <t. ^ <... be a sequence of fixed 
0 12 1 1+1 ^ 

times. Then if T is a lifetime 

p(t, ) = p{T > t. } 

1 1 

and denote the conditional probability of survival to time t^ , given 

survival to t. , by 
1+1 

P<tilt._i) = p{T > t. |t > t._^} 

p{T > t^} p(t^) 

“ p{T > ~ 
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If P(t^ = 0 f define == 0 • 



Then 




1 




(2) 



where p(t^|t^) = ' p(o) = 1 . 

Observations on Uncensored Data at Fixed Times 

Let a sample of N individuals come under observation. They are 
all observed from birth (or the appropriate event defining time zero) 
until death. V7ith the first approach, preselects a series of times, 

0 < t^ < t^ < • . • before examining the observed time of death. In the 
medical follow-up example, one might select the times corresponding to 
exactly 1,2,3,... years after a surgical procedure for cancer. An esti- 
mate of the conditional probability of survival to t^, given survival 
to t. _ is 



at time t. ,, and r. elements failed during the interval. 

1-1 1 

For a set of data which is not censored, N. = N. _ - r. _ . Now 

1 1-1 1-1 

replace probabilities by their estimates in (2) : 



1-1 




N. 

1 



(3) 



With N^ elements were present at the beginning of the interval, i.e.. 



i i 




p(t.) = lTp(tj|t._^) ( 
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N-r N-r -r 

( -) ( = — -) — 

N N-r, 



N-r -...-r. T N-r -...-r. 

_ ( 1 izi)( 1 i_) 

N-r -...-r, N-r -...-r. _ 

1 1-2 1 1-1 



= 1 - 



j=l 

N 



Now the estimate p(t^) is of the forra 



p(t^) = 



H- (r^ + r^ + . . . + r^) 
N 



and this is the same as 






N 



where is the number of the original group, of size N, that survive 

to t.. If it is assumed that the N individuals each have the survival 
1 

probability p(t), and that they die independently, then S^, the random 
number that survive to time t^ is binomially distributed, with being 
a realized value of S^. Then, considering the estimate as a random 
variable, 

;.v ■ i 



and 



N p(t.) 

E[p(t^)] = — = p(t^) 



and 



Var[p{t^)] = 



p(t^) (l-p(t^) 



Consequently p(t^) is an unbiased and consistent estimate of p(t^). 
This is true for every t^, and can be shown to be true for all t^, i=l, 
2, . . .1, as willo 
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All of this indicates that the estimate suggested is likely to be a 
good one if the sample size, N, is large. 

Clearly p(t^) • The survival probability, p(t) , is thus 

estimated at a fixed sequency of times. At each time point, t^ being a 
typical one, there are r fewer survivors than at t. ,, where r. = 0,1, 
Consequently a plot of P(t^) shows a non-decreasing step 
function, with downward steps of varying sizes at . . . . 

If the above times are close together, and if the time of death T, 
has a density function, then one can anticipate seeing values of r^ that 
are either zero or unity. 

The so-called second approach is really a limiting case of the first, 
as the time of intervals of measurement decrease indefinitely. Thus when 
a death (or loss) occurs it is only a single event. 

When no losses take place, the case now considered, the time t^ of 
the ith death is a really a realization of a random variable, denoted by 
t^; this means that p(_t^) the probability of surviving _t^, is a random 
variable. It can be shown that the expected value of p(_t^) is 



where t, < t_ < < t . 

—1 — 2 ^ 

The derivation involves integrating 



E[p(t^)] = 



CX3 

/ 



p(t) 



N! 



(i-1) ! (N-i+1) ! 



(i-p(t)]'- ^(- 



dt 



N-i+1 

N+1 



by transformation from p(t) to x; see Cramer, Mathematical Methods of 
Statistics , H. Cramer, Princeton University Press, 1946. 
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Thus one is led to use 




( 4 ) 



as an estimate of the value of p(t^), t^ being the i th time of death. 
Expression (4) provides estimator of the survival function at times of 
observed deaths when there are no losses because of censoring. The 



straight lines, or a step function with step sizes 1/(N+1) may be used. 

The estimators of equation (4) give intuitively acceptable results. 
For example, if the sample consists of only a single individual (N=l) , 
then death is equally likely to occur before or after the time at which 
the true (but unknov/n) survival function equals one half. Thus, the 
result of equation (4) is reasonable: 



The point estimates of the second approach always occur at the times of 
discontinuity forestimates from the first approach. For example, con- 
sider a data base (N=4) with deaths observed at times 1,3,4 and 7. The 
first approach gives the following step function estimate of the survival 



estimator at the points t.: t. < t^ < 
^ 1 1 2 



< t„, can be connected by 
N 



Ei;(t^)i = i 



function; 



p(t) 




0.25 



0.5 



0.75 



1-0 



0.0 



t < 7 



0 < t < 1 



4 < t < 7 



3 < t < 4 



1 < t < 3 
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The second approach gives the following point estimates 
p(o) = 1.0 
p(l) = 0.8 
p(3) = 0.6 
p(4) = 0.4 
p(7) = 0.2 

A graphic comparison of the results of the two approaches is given 
below: 




It is difficult to decide how to smooth out the step functions 
that result from the first approach. By connecting the tops of the 
"stairsteps /'one places an upper bound on reasonable estimates. By 
connecting the bottom corners of the stairsteps, one places a lower 
bound on reasonable estimates. One might draw a smooth, decreasing 
curve that passes through all (or almost all) of the vertical faces 
of the step-function estimate. The second approach suggests method 
of selecting a unique point on each of these vertical segments. 










m 



' 




Incomplete observations 



When some of the observations are incomplete, equation 
(4) requires modification. The expected value of the survival function 
at the time of the first observed death may be written: 



E[p(t^) ] 



N, 



N ^+1 



(5) 



Here is tlie effective size of the sample during the interval termi- 
nated by the time observed for the first death (o,t^^) . In the special 
case of no censoring events, the value of is unambiguous. It is 
equal to the initial sample size = N) . In this case equation (5) 
reduces to equation (4) . 

Subsequent point estimates for t^t t^,... may be calculated 
iteratively: 

1 

where t =0 and N. is the effective sample size over the time interval 
o 1 ^ 

(t. , t. ) . Thus, 

1-1 1 

i N . 

E[p(t.)] = IT (6) 

j=l j 



Variance of the estimators 

Kaplan and Meier, reference (1) , give an expression for the exact 

calculation of the variance of step functions. They also discuss 

"Greenwood's formula,” a large sample approximation that ignores terms 
2 

of order 1/N^ . 
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Herd, reference (2) , presence without derivation an expression for 



the variance of estimates using the second approach (point estimators) 



N . 



N. 



V(t ) . var {E[p(t.)l) - Tf <^7^) - TT (5^)' 

D=1 D D=1 j 



The notation here follows that for the estimating equation (6) . 
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III. THE ESTIMATORS 



This section describes the nine non-paraiaetric estimators and four 
jackknife estimators of the survival probability. It also describes 
the parametric estimator for an exponential decay function. Exponential 
life distributions are the starting point for much of reliability theory 
and practice. The estimator derived from the exponential is regarded as 
"par” when the simulated data is based on an underlying exponential decay 
distribution for deaths. Thus^ when deaths are exponentially distributed^ 
the non-par ame trie estimators may be compared relative to each other, 
and they may be compared with the parametric estimator as a standard. 

A hypothetical data base, consisting of five individuals, is used 
to illustrate each of the estimators. This sample data base is as 
follows : 



Individual 

A 

B 

C 

D 

E 



Time of Death Time of Truncation 

1 

Unknown (>2) 2 

3 

Unknown (>6) 6 

7 



The data have been arranged in time sequence of the death and trunca- 
tion events. In the medical example, the data might indicate that 
patients A, C and E were observed to die exactly 1, 3 and 7 years, re- 
spectively, after their surgery. However, B and D moved away or other- 
wise became unavailable to the observer at these times. Further, the 



cause of the unobservability is unrelated to the patient's health and 
life expectancy. 

A. STEP-FUNCTION ESTIMATORS 

1. The First Estimator , "p^(t)" 

p^(t) is a naive estimator; it is expected to perform poorly 
relative to the other estimators, p^ only depends on the data from in- 
dividuals whose deaths are observed. It ignores any information from 
the partial lifetimes noted for the censored observations. Pj^(t) is 
simply the fraction of individuals surviving to at least time t among 
those individuals whose time of death is known. It is a step function; 

t 

0-1 1.0 

1-3 0.667 

3-7 0.333 

7-co 0.00 

The naive estimator, p^(t), takes no account of the successful 
survival intervals observed for the censored individuals. Therefore 
it is biased in a downward (pessimistic) direction. 

2. The Second Estimator , ”p 2 (t)** 

P 2 (t) is the product-limit estimate® Kaplan and Meier, refer- 
ence (1) / have shown that this is the maximum likelihood estimator. The 
observed events, both deaths and truncations, are arranged in increasing 
order of occurrence: t^, t 2 f.../t^/ where N is the number of individuals 

in the sample. 

Let p(t^) denote the cumulative probability of survival of an 
individual from time zero to time t^. Let p(t|t^) denote the conditional 
probability of surviving to time t(> t^) , given that the individual has 
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already survived to time t^. Then, 






(E-1) 



If we define t^ = 0 and p(0) =1, then 



3=1 



(E-2) 



The product limit estimator is in the form of equation (E-2) with 






= 1 



If the event at t . is 



truncation 



3 






N.-l 
D 



(E-3) 



If the event at t. is 
a death ^ 



Here n^ is the number of individuals observed surviving in the interval 
t. _ < t < t.« This formulation causes the product limit estimator to 
be insensitive to the exact time of the censoring events. 

The estimator is unity from time zero to the time of the first 
event, t^, reflecting the fact tliat all individuals in our exaraple are 
observed to live until at least time t^. 

If the event at time t^ is a truncation, then the estimator 
remains at unity until at least time t^. Again, no deaths 
are observed in the sample before t^. 

If the event at time t^ is a death, then the estimator drops 
to (N-1 )/n. This drop reflects the observed death of 1/N of 
the survival sample just prior to t^. 

Values of the estimator are calculated iteratively at successive 



values of t^ (i=l, 2 , . . . ,N ) . 
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The size of the survival sample declines as truncations and 
deaths remove individuals from observation. For the hypothetical data 
base listed above, one obtains: 
t 



0-1 






5/5 = 1.0 


1-2 






4/5 = 0.8 


2-3 


(4/5) 


X 


(3/3) = 0.8 


3-6 


(4/5) 


X 


(2/3) = 0.533 


6-7 


(8/15) 


X 


(1/1) = 0.533 


7-00 


(8/15) 


X 


(0/1) = 0.0 



The product-limit estimator explicitly accounts for the sur- 
vival of these individuals (up to the time of the last death before 
each censoring event). Thus, ^ step function with a value 

that is not less than P^(^) ^or any value of t. If the sample contains 
no censoring, then and P 2 (t) are identical. 

If the last event in the sample is a truncation rather than a 
death/ then tlie modified data give the following estimate, i.e., 
individual E had disappeared from the observer at time 6.5 (so that 

the fact of E*s death at time 7 is unknown). 

^ ” Modified data 

0-1 1.0 



1-3 0.8 

3-6.5 0.533 

Since tlie time of the death for individual E is now unknown, 
one can only estimate that; 

0 ^ ^ 0.533 for t > 6.5 
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If the analyst is willing to assume a functional form for the 



survival function, then he may calculate the manner in which the 
estimator decreases to zero. However, the data alone are insuf- 

ficient when a strictly non-parametric estimator is used. 



The product-limit estimator is a useful and intuitively appeal- 



ing method of dealing with incomplete observations. It has been wider 
used and studied. However, the product-limit has one disturbing 
characteristic : 

Most of the biological, physical or other causes of deaths pro- 
duce a survival probability that continuously decreases in time. 

It is, therefore, one may be a little uncomfortable estimating 
the survival probability with a step function. One is tempted 
to smooth the estimator to make it a monotonic decreasing func- 
tion of t. 

3. The Third Estimator , "p^(t)" 



function with discrete drops at those times corresponding to the observed 
deaths in the sample population. It may also be expressed as a product 
of conditional probabilities: 



ditional probabilities on the right-hand side of Equation (E-4) differ 
somewhat from those in Equation (E-2) : 



p^(t) is a modification of p^Ct). Like p^Ct), it is a step 



i 




(E-4) 



where the t are the times of observed deaths and t is zero. The con— 
k o 




(E-5) 
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Equation (E-5) differs from Equation (E-3) in the interpreta- 



tion of the numbers of individuals at risk. Here, the value of N is 
taken to be the average number of individuals observed surviving in 
the interval between the (k-l)_^ observed death and the kth observed 
death. The number of observed survivors decrease at intermediate times 
if events are censored, and hence the are not necessary integers. 

The value of is regarded as the effective sample size for 
the interval from ^ sample data base shown cibove, 

individual B is known to have survived from time 1 to time 2 , or half 
of the interval between the first death at t=l and the second death at 
t=3. Therefore, the estimator p^ treats individual B as half a parti- 
cipant in the interval between the death of individuals A and C. 

The effective sample size for this interval is then 3.5 
( 2 - 1 ) 

(n = 3 + ^ (full contributions from individuals C, D and E, 

plus a half contribution from B) . For our hypothetical data base, the 
following values are calculated for p^t 

t 



P3(t) 



0-1 



5/5 =1.0 



1-3 


4.5 


X 


H 

0 

0 

11 

0 

• 

00 


3-7 


(2. 5/3. 5) 


X 


0.8 = 0.571 


(7) 


(1.75/2.75) 


X 


0.571 = 0.364 



The value of P 3 (t) can never be less than the corresponding value of 
p^ (t) . In the special case with no censoring events the estimators 

p^(t), P 2 (t) and P 2 (t) are identical. 

One might perturb the data by shifting the time of B's trunca- 
tion event down to l+£ or up to 3-e, £ arbitrarily small. The depend- 
ence of the estimator p^ upon the exact time of the censoring events 
may now be demonstrated. 
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For piirposes of illustration, the time of the censoring event 
for individual B{t^) is decreased from 2 to 1.1, then increased to 2.9. 



t 


P 3 (t) , = 2 


P 3 (t), = 1.1 


Pi(t), 


0-1 


0 

• 

< — 1 


1.0 


1.0 


1-3 


0.80 


0.80 


0.80 


3-7 


0.571 


0.538 


0.597 


(7) 


0.364 


0.342 


0.380 



This example demonstrates an intuitively appealing characteris- 
tic of the estimator, p^. As the total observed survival time increases 
for the individuals in our sample (with deaths held constant) , the value 
of the estimating function increases over at least a portion of its 
range. 

We may safely assume that the true survival function eventually 
tends to zero with time, since no physical or biological system lives 
forever. However, there are no observations on the survival of indi- 
viduals beyond time 7. The data only indicate that our step-function 
estimator drops to a value of .364 at t=7, but the nonparametric esti- 
mator gives no information about the survival function's subsequent 
decline from .364 to zero. However, the data alone are insufficient 
when a strictly nonparametric estimator is used. 

B. POINT ESTIMATOR 

As mentioned above, the estimators p^^, p^ and p^ are somewhat unde- 
sirable because they give step-function estimates for a continuous 
survival function. The next three estimators p^, p^ and p^ are modifi- 
cation of the first three. Again they provide estimates of the survival 
function only at those points in time that corresponde to observed deaths. 
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These estimators are specified by Equations (E-2) and (E-4) , except for 
a substitution of the term (N+1) in place of (N) o 

Since the point estimators have rigorous definitions at only dis- 
crete points in time, it is necessary to offer an interpolation rule. 

That is, we need a method of "connecting the dots." The method proposed 
here is to assume that the survival function declines in a piece-wise 
exponential decay between the discrete points in time. This procedure 
is equivalent to assuming that the hazard function is essentially con- 
stant between a consecutive pair of the discrete times, but that the 
hazard varies from one time period to the next. Such an assumption is 
intuitively acceptable unless one suspects violent fluctuations in the 
hazard function. 

1. The Estimator , *'p^(t)" 

p^(t) is analogous to p^(t) in that only those individuals ob- 
served to die are included in the sample. These two estimators are naive 
because they suppress all data from the survival times of individuals 
terminated from observation by censoring. 

These estimates, i.e., p^(t) and p^(t), tend to ignore informa- 
tion from the more long-lived individuals in the sample, and they may 
be expected to give biased estimates of the survival function. 

The point estimator P^(t) gives the following values with sample 
data base presented earlier in this section. 







Interpolation 


t 


ft 


t 


1^4 (t) 


0 


o 

• 

1 — 1 


0-1 


f £n(3/4) 
e 


1 

3 


3/4 = 0.75 
(j) X 0.75 = 0.5 


1-3 


-iln(2/3) 

(^)e 


7 


(•j) X 0.5 =0.25 


3-7 


t^.,n(l/2) 

(2)e 
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p 
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... : . .»* .. #-»£■- - ■ ' • 




f' 












4 » 





The interpolation for connecting the dots are as follows: 




t 


0 £ t < 1 


1 £ t £ 3 


3 < t < 7 


p ( t) 


e-t/T 


t-1 

p(t-l) e ^ 


t-3 

p(t-l) e ^ 


interpo- 

lation 


t = 1 

1 

3 " T 

1= ® 


t = 3 

2 

1 3 T 

1=4® 


t = 7 

4 

1 1 ” T 

1 = 1® 




- = -£n(j) 
T 4 


^ -£n(|) 


1 _ 




T “ 2 


T 4 


p(t) 


CD 

rt 

• 


t-1 0 

,3, 2 ‘‘"'s’ 

(j)e 


t-3 9 

1 4 

(j)e 



2. The Estimator, ”p^(t)" 

■ D 

The estimator p_ (t) similarly corresponds to the product-limit 
estimator p^Ct). These tv/o estimators use information from the indi- 
viduals on whom there are censored observations, p^, like does 
not exploit information about that portion of the censored observation 
after the death event (of some other individual) preceding the censor- 
ing event. 
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For our hypothetical data base the following values are 



calculated for p (t) ; 

D 







Interpolation 


t 


P5(t) 


t 


P5<t) 


0 

1 


1.0 

5/6 = 0.833 


0-1 


f £n(t/6) 
e 


3 


(|-)x 0.833 = 0.625 


1-3 


t-1 p .3 

,5, 2 


7 


(j)x 0.625 = 0.312 


3-7 


t-3 . 1 

5 4 *■''*2* 



l^enever censored observations are present, the estimator p^(t) never 
exceeds ( t) . 

For p^(t), the value of is taken to be the number of surviv- 
ing individuals in the sample just before the observation of the i th 
death. This value is smaller than the number of surviving individuals 
just after the (i-1)^ death if any truncation events occur in the 
interval. In fact, is the smallest number of surviving individuals 
observed at any time during the interval t^) . Thus p^ might be 

expected to introduce a bias by using values of that are, on the 
average, too small. However, this bias would be much less severe than 
the bias anticipated for the estimator p^(t). 

The estimators and p^ are insensitive to the precise times 
of the censoring events . A change in the time of the censoring event 
for individual B to 1+e to 3-e, e arbitrarily small, does not alter the 
estimates frora p^ and p^ given in the preceding paragraph. 
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(t)" 



3. The Estimator, "p^ 

6 

The estimator P^(t) corresponds to by accounting for all 

of the survival time for the truncated observations. For our hypo- 
thetical data base, the following values are calculated for p^ (t) : 

6 









Interpolation 


t 


P6(t) 


t 


P6(t) 


0 


1.0 


1 — 1 

1 

O 


f 2,n(5/6) 
e 


1 


5/6 = 0.833 




^ In (|^) 


3 


X 0.833 = 0.648 

1.75 


1-3 


2 '4.5 

(^)e 

t-3 n 

4 <2.75> 


7 


^ 0-648 = 0.412 


3-7 


(0.648) e 



The estimator p^ (t) is based on the average number of surviving 
6 

individuals noted in the various time intervals. These estimators give 

part credit for individuals whose lifetime is censored in mid-interval. 

The value of II. for p^ (t) is an unweighted time average. If the obser- 
1 o 

vation of an individual is truncated after 23% of the interval has 
elapsed, then that individual contributes a value of 0.23 to N^. Indivi- 
duals who are observed to survive the entire interval, and the individual 
whose death terminates the interval each contribute a value of 1.0 to N^. 
This interpretation of the effective sample size is approximate if the 
hazard is approximately constant over the interval. If the hazard function 
changes markedly within a time interval containing censored events, then 
this interpretation of the effective sample size is biased. Therefore, 
the procedure of determining the value of for the estimator p^ (t) is 
based on the implicit assumption that the survival function is locally 
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exponentialo If the hazard function may be assumed to vary slowly over 
each of the time intervals (t , t.) then p. would appear to be biased 
on an acceptable approximation. 

The estimator p^, like p^, depends on the precise times of all 
deaths and censoring events. 



Pe (t ) r = 2 






Pe(t), 



= 2.9 



0 



1.0 



1.0 



1.0 



1 



5/6 = 0.833 



5/6 = 0.833 



5/6 = 0.833 



3 X 0.833 = 0.648. (-|^) x 0.833 

7 (4m>x 0.648 = 0.412. x 0.628 



3 95 

0.628. ( 4 ^) X 0.833 = 0.665 

0.399. X 0.665 = 0.423 

2 % ! ^ 



This illustrates that an increase (or decrease) in the total observed 
survival time causes an increase (or decrease) in the estimate p- over 

D 

at least some of its time range. 

If the last event is a censored, and not an observed, death, 
these estimators also require definition for the time period starting 
with the time of the last death and ending with the time of the final 
censoring event. 

The method proposed here for p^(t) and P^(t) is to continue the 
exponential function used in the interval terminated by the time of the 
last death. This procedure can be illustrated with the modified data 
base used above in the discussion of p^ and p^. 



C. THE BAYESIAN ESTIMATORS 

Consideration is next given to quasi-Bayesian estimators based on 

a uniform prior distribution on the unit interval. Let X ,...,X be the 

1 N 

true survival times of N individuals which are censored on the right by 

N follow-up times Y ,...,Y.,. It is assumed that the X. are independent, 

1 i'J 1 
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identically distributed random variables with common distribution p(t) 
and we v;ish to estimate the suirvival function 

p(t) = Pr(x > t) 



However, we only have available the data, 



5 . 

1 



min {X^, Y^} 

’ 1 if X. < Y. 

1—1 

0 if X. > Y. , i=l 
11 






If 5 . =0, then Z. is called "a loss", and if 
1 1 

6. =1, then Z. is called "a death", 

1 1 

Then [5^_ = 1] = ^ = P(t), i=l,...,N, 

The maximum likelihood estimator for p{t) is 

N 

p(t) = — where s = Z 6. 

i=l " 

is the niomber of successful tests, s has the binomial distribution. 
P(s|p) = (^) s=0,l,...,N, 0 < p < 1 

f (p) = 1, 0 < p < 1 
The joint density of s and p is 



f 

s,p 



(s,p) 



(^) p^(l-p)“~^, 0 < p < 1, s=0,l. 



N. 



The marginal for s is 



PgCs) 



/6 p"{l-p)"-"dp = (N). 

S ^ 



s I (N -s) 1 
(N.+ l) I 



1 

N + 1 
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for s=0,l,...N. Thus, averaging over the values of p, all of which are 
assumed to be equally likely,^ the values of s are equally likely to occur. 
The posterior for p then is 




^ (N+2) S ,tl-S 

r (s+i)r (N-s+1) ^ ^ 



0 < p < 1 , 



a beta density with parameters s+1 and N-s+1. The mean of the posterior 
is (s+l)l (N+2) and the modal (maximum value) of the posterior is s/tJ ; thus 
the Bayes estimate of p (given s survivers occur in the sample of N ) is 



* _ s+1 
P “ N+2 



.-J 



(C-1) 



Then, equation (C-1) yields a step function and also has shown that the 
uniform prior has the effect of adding two individuals to the popula- 
tion at risk with one dying at time zero and the other essentially 
immortal. 

The Bayesian estimators based on a uniform prior distribution on 
the unit interval are denoted ' that correspond, 

respectively, to the estimators p^(t), p^(t). The sample 

data base thus gives the following estimates of the survival function: 



0-1 


4/5 = 0.8 






6/7 = 0.857 






6/7 = 0.857 


1-3 


3/5 = 0.6 




X 


0.857 = 0.714 


'!> 


X 


0.857 = 0.714 


3-7 


2/5 = 0.4 




X 


0.714 = 0.536 




X 


0.714 = 0.556 


(7) 


1/5 = 0.2 


(i) 


X 


0.536 = 0.268 


1.75 

4.75 


)X 


0.556 = 0.354 
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At the time of the final event (whether a death or a truncation) , 
these step- function estimators drop to some positive value. Again, 
we have no data to indicate how the survival function proceeds to zero 
at subsequent times « 

D. THE JACKKNIFE ESTIMATOR 

We will assume that we observed, or have generated in a simulation, 
a survival probagility p(tj), j=l,...,n, from various sample sizes. 
Furthermore we have some parameter or characteristic p(t^) of the 
sample size which we wish to estimate with an estimator p(t^). The 
jackknife estimator p(t,n) described below is an approximately unbiased 
estimator of p(tj). A modification of it has other useful properties. 

p ^(t,n-l) is the estimator from the sample of n of the X^*s with the 
i th value deleted from the sample. 

p^(t,n) = n p(t,n) - (n-1) p^^(t,n-l) i=l,...,n 

1 n ^ ^ n-1 ^ 

p(t,n) = — Z p. (t,n) = n p(t,n) E p (t,n-l) 

n . - 1 n . _ -1 

1=1 1=1 

the p^(t,n), called the PSEUDO-values . 

The PSEUDO-values can be used to obtain variance estimates of p(t,n) 
and to set approximate confidence limits, using Student's t. 

The idea is that the PSEUDO-values will be approximately indepen- 
dently and noirmally distributed. The jackknife estimator p(t,n) is a 

2 

sample average so we form an estimate S” of its variance given by 

p(t,n) 

the following relationship (Miller, 1974) : 

^ Sp.^(t,n) - -(Zp. (t,n))^ 

^2 1 n 1 

S = 

n-1 

si 

p(t,n) n 



30 



This procedure is particularly useful if the number of data points 
is small, but it must be used with care. Note, that the estimator p(t,n) 
is designed to eliminate a ^ bias terra in the estimator p(t,n). Of 
course the computational aspects of the complete jackknife can be quite 
onerous, especially if p(n) were, say, a complicated maximum likelihood 
estimator. liiller, reference (4) has shown that the product limit 
estimator is its own jackknife. 

Logistic Transformation 

Although one can legitimately jackknife the Kaplan-Meier estimate 
directly, there is some reason to believe that a preliminary transforma- 
tion will give improved results. Consequently, consider the transforma- 
tion 



Z = Jin (■ 



-El t ) 

l-p(t) 



•) 



and notice that where the range of p(t) is from zero to unity, the above 
transformation makes the range of % run from to o°. The procedure 
utilized will be as follows. 

(A) Compute the overall estimate at a time point t, using all N data 
points, and using a "continuity" correction that has the effect 
of removing the effect of a zero in the logarithm (see D.R. Cox, 
Analysis of Binary Data, Methuen Monograph) : 



I = £n( 
N 






-) 



(B) Compute the Ji-values by leaving out each data point in turn when 
computing p(t): for i=l,2,...,N. 



Z . . = Zn{ 

N-1,1 






* 2(N-1) 
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(C) Form the pseudo-values 



z . 

1 



n£. 



N 



(N-1) £ 



N-l,-i 



(D) 

(E) 



(F) 



2 

Compute z, 

Put approximate confidence (l-a)»100% limits on E[£] 
L <_ E[£] < H 

where H (L) = z +(-) t^_^ (N-1) / ^ 

Transform bash to obtain 

L H 



, and 



1+e 



1+e 



H 



as follows 



The true value, p(t), should be enclosed between these levels for 
roughly (l-a)*100% of all samples. The coverage properties of this pro- 
cedure will now be checked by simulation; successive samples of size N 

will be selected, the jackknife limits H and L will be computed for each, 

L H 

e e 

and a check will be made as to whether — < p(t) < — or noto Tables 

1+e^ 1+e^ 

illustrating performance are given subsequently. 
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Let 



N' 



is the logistic transformation estimator from the 



of the X^*s with the ith value deleted from the sample. 



Vi,-i ■ 






1-p , . (t)+— i 






i 

t 


1 


2 


3 


4 


5 


h 


3.04 


0.98 


0.98 


0.98 


0.98 


^2 


3.04 


0.98 


0.98 


0.98 


0.98 




0.63 


0 


0.98 


-0.46 


-0.46 




0.63 


0 


0.98 


-0.46 


-0.46 






-3.04 


-3,04 


-3.04 


-1.89 



z- = - (N-1) , . 

1 N N-1,-1 



= N£n( 



+ 2N 



2N 



r-) - (N-1) 3 






■) 



1-p, 



(t) + 



N-l,-i' ' 2 (N-1) 



(ll) are called PSEUDO- values of logistic transformation j 



following values are calculated: 



i 

t 


1 


2 


3 


4 


5 


h 


-6.05 


2.198 


2.198 


2.198 


2.198 


^2 


-6.05 


2.198 


2.198 


2.198 


2.198 




-1.9 


0.606 


-3.314 


2.446 


2.446 


^4 


-1.9 


0.606 


-3.314 


2.446 


2.446 




-3.0626 


-3.0626 


-3.0626 


-3.0626 


-7.162 



sample n 



the 
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Average of the pseudo-values 



z 




z . 
1 



Invert to find jackknife estimator of logistic transformation 

Z = iln 

1 , z -i_ 

p(t) = jjCTle .. - 2N 



called the jackknife estimator 
of logistic transformation 



Variance of the 



S 



2 



z 



Var ( z ) 



1 

n-1 



n 

Z 

i=l 



■ 



z 



The following values are calculated: 



t 


z 


p(t) 


Var J t 


h 


0.5484 


0.646 


13.6 




0.5484 


0.646 


13.6 




0.0568 


0.516 


6.727 


^4 


0.0568 


0.516 


6.727 


^5 


-3.882 


0 


3.361 



The jackknife estimator for estimating variability and giving confidence 
interval. 

Tukey, reference C3) has suggested that in the jackknife procedure 
we consider tlae pseudo values Z^(n) as approximately independent and 
identically districuted and consequently, since z is an average of 
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I H I t4i JSJS®'? 












(Miller has shown) 



the Z ^ (N> , proceed as if 



N ^ z - £ 



N 



1 ^ _ 9 1 ^ 
z (z. - z, = )- 
1=1 



has t-distribution with N-1 d.F . 

If the are approximately normal variates 
confidence bands for the unknown p(t) are given, as for the mean of 
any normal variate when estimated from sample size n. 



z + — t. (N-1) 

- 1 - 0/2 



(D-1) 



1 . e. 



2 - — 

i-«/2 



p(t) + — s 

(N-1) < £n( ^) < i + — t 



l-l5(t)+ i 






(N-1) 



L(n) = z - t 






1 , L(N) 1 1 X L(N) 1 

2N ® ” 2N ^ ^ 2N ® ~ 2N 

< p(t) < 



1 ^ L(N) 
1 -I- e— 



1 + e 



L(N) 



The following values are calculated; 





4 


t, = 2.776 

l-a/2 


Lower Int. 


Upper Int. 


h 


0 


1.0 


^2 


0 


1.0 


S 1 


0 


1.0 


^4 


0 


1.0 




0 


0.14 
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The basis for this leap of the imagination seems to be that if 
^ = X = then the procedure for obtaining confidence intervals 
using equation (D— 1) and pseudo— values is the same as the procedure 
using jackknife. Then if = z and 



2 




Z . 
1 



we have 



2 . 
1 



= N£ 

N 



= N X. 



N 



(n-1) 






I 

N-l,-i 

i''! 

{ Z X. } - X. 

j=l ^ 

N-l 



N 

= Z X 



N 

[ Z X.] + X. 

j = l ' 



X. 

1 



Thus the pseudo value 

z. = X. and z 
1 1 

The pseudo values are independent if ^ ~ they are noirmal if 

X. is normal. 

1 

E. PARAMETRIC ESTIMATOR, "p^ (t) ” 

This paper considers one additional estimator, denoted p^ (t) . It 
is a parametric estimator. Therefore, it is not really a competitor to 
the thirteen non-parametric estimators considered here. In general, a 
parametric estimator would not be used if the functional form were re- 
garded as unknown. Similarly, a non-parametric estimator would not 



1 ^ - 
= - Z X. = X 
n . _ 1 n 

1=1 
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normally be used if the survival function were strongly suspected to 
have a specified form. 

p^(t) is the well known maximum likelihood estimator for the expo- 
nential distribution : 






= e 



-t/T 



where 



E t. 

1 

number of observed death 



In our sample data base, the total observed survival time is 19, and 
three deaths are observed. Thus, 



and 



T. t. 
1 



= 1 + 2 + 3 + 6 + 7 = 19 



T = 



3 



P^(t) 



= e 



-3t/19 



Calculations for selected times of interest yield the following esti- 
mates : 



p^(0) = 1.0 

p^(l) = 0.854 
p^(3) = 0.623 
p^(7) = 0.331 

The thirteen non-parametric estimators are compared for a variety of 
generating distributions for both the death mechanism and censoring 
mechanism. 
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IV. INSTRUCTIONS FOR USING PROGRAM 



INPUT 

Each input card bears nine variables. The distribution of time of 
death is entered in the first set of (five) columns, the censoring dis- 
tribution is entered in the second set of (ten) columns, a parameter of 
the censoring distribution is entered in the third set of (ten) columns, 
the number of replication is entered in the fourth set of (five) columns, 
the niomber of the event is entered in the fifth set of (five) columns. 
For the purpose of all print output used code "0” and *'l*' in the sixth 
set of (five) columns, the seed number is entered in the seventh set of 
(five) columns, after the card giving the time of the last event of a 
data set, a card with **0" or "1" in the column 50 is inserted, i.e., the 
"O'* indicating more data sets to follow and "1" indicating the last data 
sets and t value is entered in the ninth set of (eight) columns. 

The distribution of timeof death and of censoring time used code as 
follows; 

Code Type of Distribution 

1 Uniform 

2 Exponential 

3 Delta function 

OUTPUT 

The output lists: 

1) the time of each observed failure 

2) estimated survival probability at that time 

3) the variance of that estimator 

4) result of goodness fit 
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a) mean error 

b) mean absolute error (ABS) 

c) root^ me an- square error (RMS) 

5) total number of observed death 

6) confidence interval at particular time 
Definition of Fortran Variables 



NDIE 

NTRUNC 

XTRUNC 

NREPL 

NEVENT 

NWRITE 

NEND 

TN 



10 



11 



'12 



the distribution of time of death 
the distribution of censoring time 

the parameter of the distribution of censoring time 
number of replication 
number of event 

write all output or partial output of simulation 

indicate more data sets or last data set 

t statistic value 

the estimator, p^(t) 

the estimator, p^ (t) 

the estimator, 

the estimator, P^(t) 

the estimator, p^(t) 

the estimator, p (t) 

6 

parametric estimator , p^(t) 
jaclcknife estimator of logistic transformation of P^(t) 

jackknife estimator of logistic transformation of P^ (t) 

jackknife estimator of logistic transformation of p (t) 

D 



Bayesian estimator of p^^Ct) 
Bayesian estimator of p^ (t) 
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Bayesian estimator of 
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^14 


: jackknife estimator of logistic transformation of (t) 


SL(I,J) 


: PSEUDO-value 


SBAR 


: average of pseudo-value 


Var 


: variance of estimator, p(t) 


Var_ 

J 


variance of jackknife estimator 


u(I,J) 


mean of goodness fit 


w(I,J) 


absolute mean of goodness fit 


s (I,J) 


root mean square error 




: upper confidence interval of P ^_4 


s 


: lower confidence interval of p, , (t) 

14 


s 


: upper confidence interval of P 0 (t) 


^4 


: lower confidence interval of p (t) 

o 




: upper confidence interval of p (t) 




lower confidence interval of p^ (t) 


S 


: upper confidence interval of 




: lower confidence interval of 
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To compare RMS with product limit (p (t)) and jackknife estimator of 
logistic transformation (p^^(t)) 



1 


A.OOOC ICCO 


20 1 


505 : 


Input 










1 


C.C127S 


C.95C00 


C. 93 726 


• 


Output 










c 


0.01581 


0.0 


C.COOOO 












1 


C. C3C8S 


C.S9722 


C. 88376 














1 


C.CS^3C 


C. 8^^^^ 


0.83199 














1 


C .11S25 


C. 79167 


C. 78256 














1 


C . 1 i £ 7 7 


c. 73€es 


C. 73445 














0 


C. 17t7^ 


G.C 


0.05254 














c 


C. 1ES82 


c.o 


C. 00000 














1 


C.ISC^^ 


C. 6 7 731 


0.68724 














1 


0.1S81 6 


C. 61574 


C. 63201 














1 


C.2767C 


0.55417 


0. 57755 














1 


C.2HC60 


C. 49259 


C. 52356 














a 


C. 21870 


0.0 


0.00000 














1 


C.^EC75 


0.42222 


0.47045 














1 


0.6^596 


C. 35185 


0.40828 














1 


C.7C112 


0.28146 


0.34624 














1 


c. / /ess 


C. 21111 


0.26413 














1 


1.18657 


0. 14C74 


0.22201 














0 


1.25112 


0.0 


0.05254 














, 1 


1.^7370 


C.O 


0.05254 














C.IOC 


C.20C 


0.300 


0.400 


0. 500 


C.600 


C.7C0 


0.800 


C.900 




-0.0C2 


O.OOC 


0 .002 


-0.000 


-0.001 


-0.001 


0.007 


0.028 


0.066 


N6AN 


-0.013 


-C.007 


0.C04 


0.015 


0. 029 


0.044 


0.069 


0.109 


0 . 166 


WEAN 


C.053 


C.C71 


0.C84 


0*090 


0.097 


0.093 


0.091 


0.077 


0.078 


ABS 


0.057 


G.G68 


0.075 


0.080 


0.067 


0.087 


0.094 


0.116 


0.166 


AOS 


0.06^ 


C.091 


0.105 


C.l 13 


0. 119 


0.115 


C.113 


0.100 


0.102 


PWS 


0.071 


0.085 


0.C94 


0.101 


0. 108 


0.109 


0.120 


0.140 


0.185 


RWS 


C-965 


C.951 


C.901 


0.841 


0.781 


0.722 


0.661 


0.678 


0.751 


CCNF 


C.767 


C. 591 


0.459 


0.355 


0.267 


0.185 


0.116 


0.057 


0.005 


CCNF 


57.726 


9 1 .25^ 


97.065 


97.959 


93. 542 


99.125 


97.959 


98.542 


97.668 


PER 


lOGC 


lOCO 


ICCO 


lOOC 


998 


990 


951 


796 


343 





2 


2 


<I.OOOC loco 


10 1 


1509 : 


Input 










1 

1 

1 

1 

1 

1 

1 

I 

1 

1 


C 

1 

0 

1 

0 

1 

1 

1 

1 

1 


C. 03688 
0.14610 
C. 22306 
0.22401 
0 .30^47 
C.4635S 
C. 71228 
0. /950C 
1.1^699 
2,27255 
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Computer output of the fourteen estimators 
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V. RESULTS OF THE SIMULATION 



This paragraph presents graphical comparisons of the 13 estimators 
based on simulated data. Each comparison is based on 1000 replications 
of a simulated data base. The bias and RMS error (square-root of mean- 
squared error) of each estimator depends on the parameters that control 
the simulated data base. No single estimator dominates all others 
under all conditions. 

The bias and RMS errors of the estimators depend on several factors: 

(A) The sample size (NEVENT) of individuals under observation at time 
zero affects the accuracy of the estimators. In general^ a larger sample 
size leads to a better estimate than a smaller sample. Values of NEVENT 
selected for simulation are 5, 10, 25, and 50 (plus one simulation with 
NEVENT = 100) . 

(B) The distribution of times at which the observations are censored 
(unless the individual dies earlier) affects the performance of the 
various estimators o This distribution is particularly important in con- 
junction with the distribution of lifetimes (do most individuals die 
before censoring is likely?, are deaths and censoring events about 
equally likely at all times?, are most observations censored before death?) . 
Three types of distributions are assumed to underlie the censoring mech- 
anism: 

(1) Some of the samples are generated on the assumption that no censor- 
ing occurs. 

(2) Some samples are generated from a uniform distribution of times of 
censoring. 
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(3) Other data bases are generated from an exponential distribution of 
censoring times. 

(C) The distribution of lifetimes (ignoring the possibility of censor- 
ing) also affects the performance of the various estimators. Two types 
of distributions are assumed to underlie the death mechanism: 

(1) Some of the samples are generated from a uniform distribution of 
lifetimes. 

(2) Other data bases are generated from an exponential distribution of 
lifetimes . 

If a uniform distribution of lifetimes is selected, its range is 
always over the interval from time 0 to time 1. If an exponential dis- 
tribution is selected, it always has a mean lifetime of 1. The distri- 
butions of truncation times (uniform or exponential) have parameters 
.25, .5, .667, .75, 1, 1.333, 1.5, 2 and 4. A wide variety of samples 
may be simulated by mxing various pairs of distributions (for censoring 
times and deaths) . Since the time units are arbitrary, the restriction 
on mean lifetimes is irrelevant. 

The true value of the survival function is, p(t), and the form of 
this function affects the relative performance of the 13 nonparametric 
estimators. For example, the Bayesian estimator tends to be 

better as measured by square-root of mean- squared error than its counter- 
part (the product-limit estimator, p^(t)) for the time frame in which 

.3 < p(t) < .9 

However, the product limit estimator tends to be better for those times 
when p(t) is close to zero or unity. 

The point estimators, P^(t) and P^(t) tend to be better than the 
product-limit estimator (p^(t)) for all time periods. The jackknife 
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estimators of logistic transformation (Pg(t), p^(t) , Pj^g(t)) of point 
estimators tends as same as its counterpart point estimators (p^(t), 

Pg(t), p^(t)) for all time periods. And also the estimator fomed by 
jackknifing the logistic transformation (p^^(t)) of the product limit 
estimator tends to be better than its counterpart product limit (p^ (t) ) 
for the time frame in which 

.1 < p(t) < .7 

However, the product limit estimator tends to be better for those times 

when p(t) is close to unity. Point estimators, p- (t) and p^ (t) tend to 

b 6 

be same for the time frame in which 

0.1 < p^(t)< 0.9 

However, the Pt-(t) tends to be better for those times when p (t) is 
close to unity. The jackknife procedure may be validated, in an empiri- 
cal sense, by sampling experiments or computer simulation in the follow- 
ign manner. First, times of censoring and death are obtained by drawing 
random numbers from postulated distributions. Second, the jackknifed 
estimator of the logistic-transformed product-limit estimation is found, 
and confidence limits are computed by the method of Tukey, reference (3) . 
Since the true value of survival function, p(t) , is known, so is the 
theoretical value of A. The jackknife confidence intervals can be 
checked for coverage: if _< A ^ then the particular interval covers, 

while otherwise (if A < or < A) it does not cover. Finally, the 
above procedure can be repeated many times (say 1000) and the fraction of 
repetitions which contains the true value of A is recorded. This fraction 
of the coverage should desirably be close to (1-Ot) , the nominal confidence 
level. 
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The jackknife confidence limit procedure can be said to be robust 
of validity / ref (7), if the actual coverage is close to the nominal 
coverage, 1-a, for a various distributions. Such seems to be true for 
large n (n ^ 50) . However, the jackknife confidence limits do not cover 
accurately when the true value of p(t) is close to unity. 

The following tables illustrate confidence limits of jackknife 
method of product limit (p^Ct)). Many computer generated graphics are 
presented on the following pages to complete this section. 
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Simulation Experiments Validating Table 
95% Confidence Limits (t value = 2.776) 
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simulation Experiments Validating Table 
=10 95% Confidence Limits (t value = 2.262) 
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simulation Experiments Validating Table 
=25 95% Confidence Limits (t value = 2.093) 
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simulation Experiments Validating Table 
95% Confidence Limits (t value = 2.010) 



o 

in 







50 



NO : 

I • 

. N 3 '»L “ 0 1 



L 

> • V s . * 



: ;i .0 . = : ;c 

: ( i- ; . .' s . : I ^ ) # .vs. t’ ( y ) 



Nr v» \ T- 






; \ A ^ 



« 

<> 

« 

* 






, .'rOl L-Cl 



2. AClE-01 



2.1 t<*t* )l 



2 . 



l.r.B6n-01 



1 .BS6E-CI 



1 .t d9L-01 



1.62SE-C1 



I .3 lU -01 



1. 371L-01 



I .lUF- n 



V »? f « V 

J. ‘J 



^ ^ V ? ; C ♦ >-• =* V ••• ^ V •. *f ■• 






-v< * + ♦ 



- if.*. « <• 4 - > 



V.) • 



J. 



+ 

s 4-’? *5^ 

.CO 



1. llrC-Ol 



Fig. 1: Comparison of RMS derived from sample size 5. 

(Step function estimators vs. point estimators) 



51 



L - 2 /} A -- 



'r . . C C r •• V " T - 



V z { • i H ^ ^ j » / > • ( -i ] • V s , r* I c { *•' j 



o,j 0.:^ 

««.»<•: / > ' < « t: f it i ♦ J* ^ 

Z,lLJ^~Ji *■ 



(’.‘-0 o./d 1. 

^ i '. t r s: < f /• » ^ If, «r ^ <;« 4 v x ^. < 4- «r « 



^•0 

ft « 



, uci“ 0 l 



i.^7/':--)i f 



1 . C77P-CI 






1.7C3F-;iJ 



1.610*^- Ji 



■ 1 CF-01 



1. ♦27t~ijl 



I . ATTr-C! 



I .2^?5- J1 



l.?^3L-01 



l.*}oJ?-Oi 



=> ; i- ^ 1; Jt - 

J.O 



f ^ if :« »- 4: ;-• f X X 



C-. V V t >!' < 4 4« »;« « t tf' <•• ❖ X Vf h X« ^ f < : 

0.?*^ O.IO 



0. 7S 



1. no 



1 . C7 cr-ci 



Fig. 2: Comparison of root mean square derived from sample size 5. 

(Step function estimators vs. jackknife estimators of 
logistic transformation) 



52 



4 



m;m- I 



PT ( . ) . VS . p: t 



A I '.Of. 



Ni. vrf.n - 



> I .vs. I> I ; ( ' ) .vs. 1 3 ( < ) 



I • o 

If c v>r ^ sc * 



r>. jO 






l.OO 









-v)J 



2 . 



r^oir-oi 



♦ ?.ACir-Ol 
« 

« 

« 






♦ 

« 

<f 

« 

« 

« 

« 



2. l^AE-Ol 



.3 3«.t-Ol + 

If 

if 



1. Ff /’.*'"Cl 



L.62VF-01 



l.-WK-Ol 



l.^;2SF-01 



1.27lf-01 



I .1 L3r-Ol 



# <. 4 C t ± ^ ' ! 

0.0 



; 4 < < < < - - 



'• < - 1 ; 

c. 



t 4 i. C 4 <s <- ❖ V 






n.so 



0.: 



I « 4- Is C 

l.OO 



1. 112F-01 



Fig. 3: Comparison of root mean square derived from sample sizes. 

(Step function estimators vs. Bayesian estimators) 



53 



N L> ; I = i M n I t 



.jf.i 



A . u; J«’ NT vr rj T n 



pj i . ) .v:>. .\s. !’:(■.), ,vr>. I’^O'j 



1 .6i>‘3r-oi 



5 ^ ♦ ♦ ^ ^ V ? «! V 



0 . ‘ 
' \ A X > t A 



{. < v‘. < < ■> <i * 1 



O . > 0 \j » I j 

A » <<>«<: t V ♦>!« f « « t ^ ^'s , <^ <£ ♦ J: f 4' t < < i 



I .00 

K:* 

+ 



l.^ 



i.»j32F-01 



1. 532f-Cl 



I . 3 vSC - ul 



1. i SSP-Oi 



l.’60C-Ul 



1. 26CC-01 



1. L33e-01 



1. 133‘=-Cl 






<:.9^7r-02 



i . 3 o 6 r - j r 



^ < ♦ sM < 

0.0 



< 4 V X ♦• V * 



*t If i. if * 4 4rf*A/ ■ 

o.?*> 



❖ 

+ 

0.3U 0.75 l.OO 



f .66<e-02 






Fig. 4: 



Comparison of root mean square derived from sample size 10. 
(Step function' estimators vs. point estimators) 



54 



M I hi-r - 



. . I r 



10 



p;m . ) .vs. i» ^ { «■ ) . V '» . ) .vs. p I 0 ( / ) 



.0 

.-I- N V .11 V- ? 



I .^7Ut-0l 



■ f ;i: 



J.' ‘ 



0 . ' 0 

: * < < < f V f r- < <: < f » < 

i 



0./') 1, .f;0 

VV f < < C ■^ .'• >!« «m: ^ < 1. K- f*' ^ K> ^ ^ n< ^ 



1.6701-01 



I . •? 



n 



l.^40r-oi 



\ lOL - Jl 



1.410F-01 



1.23JF-01 



1.23Cr-Ol 



1 . 1 5u - -01 



1. \^C^-C) 



I . J2oF- Jl 



i.c-^ri -oj 



1 .VOvTt-- 



V V t ^ ^ 



i . ) 



^ 4 ^ / |~ 4 

0. ,'S 



: 4 s; f 4 C- 4 ^ 

o.so 



;- jr.- V* 4: f t i 4 »? « t- f <•' /' * ap < < < ?: « + C t V'^ « < * ' 

0.7b l.OO 



e.^-ccf -0? 



Fig. 5: Comparison of root mean square derived from sample size 10. 

(Step function estimators vs. jackknife estimators of 
logistic transformation) 



55 



I 



I 1 If I - J 






^ . C C c. 



.v^ Vt N1 ' 



10 



I • * 



r’ . I ♦ ) . '. ^ • 



. / \ • ) 



- S , IM 3 ( a ) 



i.oose-oi 



J .0 



> f r 



O.UO 

^ <4 Ur 4. - > ■ <: ♦ V X X <■ »■;'«;»» X i?. X < >;•• ^ ♦ « X. # UtUc U : 



). 7S 

^ X 4 X; X <'/=>< <; » <■ 4- <■ > 



1 . 00 
: i - A ♦ t <c 






l.e>32t-OI 



I. 









U26^t- Jl 



i.refE- 



l, 133E-01 



K1 ??f-Ol 



1, ? >?r- u 






1 , of > 6 L -02 



0,0 



, X: t 4 X: X X- *!' A X- X < X- fX- * X X V X 

0 , 75 



u ;><: 4 -<; -t 

1 . CO 



f i f^-Zl 



Fig. 6: 



Comparison of RMS derived from sample size 10. 
(Step function estimators vs. Bayesian estimators) 



56 



1 4T..Lfc- ? A.occr^ N»vr^‘. 

P2(.) .vs. HS(M ,vs. .vs. >>MX» 



J . ) 

t-i 






v). ' 



I .i) t<»K-Jl 



r ;, ♦ < - r- 1 ;* 1 1 ^ t <; i;j •; t « f * is t <• 4 c < ^ i- « ^ 



0. 7b 

*•<«< <-s <>>4 /. <r:;i>;iV<»<i<sf4't < <4 i< 



I .( ^ 

i < J> 



I.CS4F-01 



>. Sj 21 - 3 £ 



S- 6t 2F-02 



». -I 



£. FSSF-C2 



J. I 15F.- D? 



6. U5F-02 



^.i^lE-D2 



7.3^ir-J2 



-» . 1 J2 



C .56&F-n2 



>. 7 7^»C- J2 



4 V > 

0.0 



0.2S 



t - 

0.‘»r C.75 l.Cu 



5. 7C4iF-02 



Fig. 7 1 Comparison of RMS derived from sample size 25. 

(Step function estimators vs. point estimators) 



57 



I 



= I 



'■» - 1 > rsliv: i: I ^ 



Pi, 



U >*•>.- -v; I 



J . »’. ( 1 - J . / ,, , i ^ ^ 

■^ • -i 0 *' • 

<• t i }. ;. c: ? ? > (! f j ^ ^ ^ ^ ^ ^ ^ ^ ^ 0 , ■ > ( I 

r ' '” <<-vvv5^;; ’<t*, 



I J .VS, '■' I c ( V , 
0. '‘3 






I . '<o 

« ' ’'r < 4 ❖ t 



1 . “ 5CI -Cl 



1. I t7 - ,1 



1 . I 7‘'-Oi 



I .04 -01 



? .400^- J2 



.367«-- )2 



333r-02 



l.C 4 ?r-oi 



*? • 4 cc r -0 / 



•3. ?( 7F-02 



7.33:^-c^ 



3 0wf - ■): 



- <- + ^ V 
) . 0 



^ ^ 4 / < < I V < < - < ^: V .1 



J:** /: 



V4 ^ Ki if t. ^ y_ 

o.so 



/ l.rjo 



f . . 30Cf.-0 3 



Fig. 8: Comparison of RMS derived from sample size 25. 

(Step fiinction estimators vs. jackknife estimators of 
logistic transformation) 



58 



L= I n f.. 



V I = 






p ? i • ) . V H ( ♦ ) . V j • »’ I ^ i v= ) . y s • p 1 1 ( a i 



^1^-01 



<)• >) 

t. <r i «.> ■:»!! 






. 1 . r S 

< A 4 4 >; ❖ . 



V j » 7 ^ 100 



1 . CAAt-C I 



^28- )a 



3‘J^-02 






6. '=ecE-C2 



l*J8-02 



0. U5P-C2 



AlP-02 



?• 2A1E-0? 



038-02 



6.560^-0: 



VA8-02 



^ V r t 4f Jt 

J ,0 



, 7SAP-C2 



Fig. 9: Comparison of RMS derived from sample size 25. 

(•Step function estimators vs. Bayesian estimators) 



59 



APPENDIX A 



ESTIMATORS FOR GROUPED DATA 



With grouped censored data the definition of p(t^/t^ given 
by equation (5) does not hold unless the assumption is made that all 
truncations occur at the end of the time interval. If, on the other 
hand, it is assumed that all truncations occur at the beginning of 
At^ the equivalent form of equation (5) is 

N. - a. - r. 



P(t^/t 



i-1 



) 




(G-1) 



Witli N^ elements were presents at beginning of interval, i.e., at 

time t. ., r. elments failed during the interval, and a. elements 
1-1 1 ^1 

truncated from the sample during the interval but prior to failing. 

As a hypothesis, assume that all aborts occur simultaneously somewhere 
within the time interval, so that r* failures occur prior to the 
truncations and time remaining r^ - r* after the truncations. Then 



N. -r’ N. -a. -r. 

VVl^ = • n' - r' - a ' ■ 

1 111 



(G-2) 



Thus, the value of p(t^/t^ depends on when the truncations occur. 

It is assumed that this is not known for the grouped data case. Never- 
theless, it is possible to place limits on the value of p(t^/t^ 
since equation (G-2) always gives values between those of equation (5) 
and (G-1) . Thus 



N . - a . - r . N . - r . 

- P' Wi> i ^n'" " 
11 1 



(G-3) 
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For average sample size approximation, a simpler expression from the 
point of view of computational ease may be derived by substituting 
a/2 for a in equation (G-1) giving 
N - — - r 

P(t. ) = (G-4) 

U - - 
2 

The equation (G-4) may be thought of as the result of assuming that 
the average number of elements in the time interval is the number at 
the beginning decreased by half the number of truncations. 

Records are usually available to provide a fairly precise time the 
deatli events. In the medical example, the exact time of death is 
usually recorded in medical records required by law. In the equipment 
lifetesting example, the time of malfunction or failure is usually 
known very precisely if the results are catastrophic; and maintenance 
records give a reasonably precise time even if the failure is not 
critical to a larger system. In the military example, the event of 
interest is usually a sensor detection or some other action that is 
routinely recorded in a log book, 

Equaiton (G~4) is a modification to the product-limit estimator, 
p^f when the times of truncation are known only in grouped form. Herd, 
reference (2), suggests a similar modification to estimators using the 
second approach (p^ or p^) \\rith aggregated truncation data. Illustrate 
results for this method based on the sample data base of the main test 
are given below. Here, of course, we do not know that individual B 
dropped out of observation at time 2 and that individual D dropped out 
at time 6. We know only that the two truncations occurred in the 
interval (1,3) and (3,7), respectively. 



61 



Product limit's modification is denoted by and Herd's modi- 

fication is denoted by p^' (t) . 

Their results on the sample data base are as follows. 



t 

0-1 

1-3 

3-7 

(7) 



P2’(t) 



5/5 = 1.0 
4/5 X 1.0 = 0.8 

2. 5/3. 5 X 0.8 = 0.571 

0.5/1. 5 X 0.571 = 0.190 



0 

1 

3 

7 



P3' (t) 

1.0 

5/6 X 1.0 = 0.833 

3. 5/4. 5 X 0.833 = 0.648 
1.5/2. 5 X 0.648 = 0.389 
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APPENDIX B 



LISTING OF COMPUTER PROGRAM 



C 

C 



C 

C 

C 



C 

C 

C 



c 

c 

c 




1 FCRMAT (15, 
i FCRMAK IX, • 
i FCRMAKiX,' 
A FCRMAK IX, ' 

5 FCRMAKIX,' 

6 FCRMAKIX,' 

7 FCRMAT (3X, 
e FCRMAT (IX, 
9 FCRMAT (IXi 

1C FCRMAT( ' • ) 

11 FCRMAT( ' 1 ' ) 

12 FCRMAT( IX, I 

13 FCRMAT (IX, 

14 FORMAT (IX, 

15 FORMAT (6X, 

16 FCRMAT (IX, 
18 FCRMAT (IX, 



I10,F10.4,5I5,F3.3) 
NDIE ERROR') 

^TRUNC ERROR' ) 

NREPL ERROR' ) 

KEVENT ERROR' ) 

NWRITE 6RRCR' ) 

9F8.3) 

2I5,14F8.5) 

'P' ,I 2,9F8.3,3X, 'MEAN 



5,10F1C.5) 

2I5,F10.4,4I5 ) 



•P' , I 2 ,9F8 
9(14, 4X) ) 
•P' ,I2,9F8 
215) 



,3 ,3X, ' ABS ' ) 



,3,3X, 'RMS') 



READ INPUTS AND SET INITIAL VALUES 
A=C.01 

"5 READ(5, l)NOIE,NTRUNC,XTRUNC,NREFL,NEVEiNT ,NWFITE,ISEEC, 
^NEND,TN 

l!piTElb?i3 ) NDIF , NTRUNC , XTRUNC , NR E F L , NEV ENT , NWBITE, ISEE 
^ D 

"write (6,10) 

NF = 14 

SRH=SQRT( .5 ) 

DC 26 1=1,1000 
T( I )=0. C 
Pl( I)=0.0 
P2( I)=0 .0 
P3( I )=0.0 
P4(I)=0.0 
P=(I)=0.0 
P6(I )=0.0 
TJ( I)=0.0 
Pil( I )=0.0 
P12(I )=0.0 
P13( I)=D.O 
26 PZ( I )=0.0 

TEST INPUTS 

IF (NDIE.lt. 1) GOTO 105 
IF (NDIE.GT.2) GOTO 105 
IF ( NTRUNC . LT.l ) GC TC 103 
IF (NTRUNC. GT. 3) GOTO 103 
IF (NREPL. LT.l ) GO TO 104 
IF (NREPL .GT. 1000) GOTO 104 
IF(NEVENT.LT.2) GOTO 102 
IF (NEVENT.LE.IOOO) GOTO 200 

ERRCR MESSAGES 

102 WRITE (6,5) 
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OOO OOo ooo ooo oc^o 



STCP 

103 ^vRITE ( 6t3) 

STCP 

1 C4 WRITE (6,4) 

STCP 

1C5 WRITE ( 6,2) 

ICE STCP 

START MAIN CALCULATICN 

2CC CC 250 J=l,9 
NN( J) =0 
CC 250 1=1, NP 
S (I ,J)=0 
U( I , J) =C 
C( I,J)=0,0 
25C W(I,J)=0 

DC 4999 IREFL=1 ,NREPL 
NCI=0 

DC 999 IEVENT=i ,NEVENT 

CREATE TTRUNC( ) 

CALL RAND0M( ISEED,TTR,1) 

GCTO (300,350, 400 ) tNTRUNC 
GCTO 103 

30C TTR=TTR*XTRLNC 
GCTO 500 

35C TTR=-XTRUNC*ALOG( TTR ) 

GCTO 500 

4CG TTR =XTRUNC 

CREATE TDIEO 

5CC CALL RAND0M(ISEED,TDI,1) 

GCTC (800,700) ,NDIE 
GCTO 102 

700 TCI=-ALOG(TCI) 

DETERMINE SMALLER CF TDIEO ANC TTRUNC() 

8CC IF (TOI.LE.TTR) GOTO 810 
TT=TTR 
ITT=0 
GCTO 850 
81C TT=TDI 

NCI = NDI +1 
ITT=1 

85C T(IEVENT)=TT 
IT( IEVENT)=ITT 

CRCER DATA 

67C IF ( lEVENT. EQ.l ) GCTO 999 
I I=IEVENT-1 
CC 890 1 = 1 , II 
IF (TT.GT.T(I) ) GOTO 890 
I 11=11-1 + 1 
CC 880 J=1,III 
JJ=IEVENT-J 
IT( JJ+1 )=IT(JJ) 

38C T(JJ+1)=T(JJ) 

IT( I )=I TT 
T(I ) = TT 
GCTO 999 
8SC CCNTINUE 

999 CCNTINUE 

TN=T{ NEVENT ) 

T7 = 0 

DII=NDI+2. 

I F (NDI .GT.O) GOTO 8122 
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(l\EVENTJ = SRH 
P=(NEVENT]=SRH 
TT = 0 

DC 3121 J=ltN£VENT 
8121 TT=TT+T{J) 

DN=TT/T ( NEVENT ) 

P£(NEVcNT) =SQRT(DN/( CN+1.) ) 

GCTO 1111 
C 

C CALCULATE Pl() AND P4() AND PIK) VECTORS A^C P7 DATA 

C 

8J.22 N=NDI 
J = 0 
I 1 = 0 
I 1 1=0 

DC 2199 I=ltNEVEMT 
T7=T7+T ( I ) 

IF( IT< I ).EQ.l) GQTC 2150 
J = 1 

GCTO 2199 

21 =C Pl( n=FLOAT(N-l)/FLQAT(NDI) 

P4( I )=PL0AT (N)/FL0AT(N0I+1) 

P11(I)=FLCAT(N)/DII 
TTT=T( I ) 

I II = II 
I I = I 
J=0 
N = N-1 

2199 CCNTINUE 

T?=-NDI /T7 

IF (J.EG.O) GOTO 2221 
TI I=T( I I ) 

CTI=T(I I n-TI I 
IF (NOI.GT.l) GOTC2200 
TAU=-TI I/ALCG( P4( II) ) 

TTAU=-TII/AL0G(P11(II )) 

GCTO 2210 

22CC TAL= OTI/ALCG( P4(I I ) /P4( II I ) ) 

TTAU=(T( I I IJ-TI I )/ALOG(Pil( ID/PIKIID) 

221C DT=TN-TII 

IF (DT.GT .150'TTAU) TAU = 0T/150 
I F(DT .GT.150*TTAU)TTAU=OT/150. 

Pli(NEVENT) =Pii ( I I )*EXP(-DT/TTAU ) 

P4 (NEVENT) =P4( I I )*EXP (-DT/TAU) 

2221 IF (NDI.NE. NEVENT ) GO TO 2225 
DC 2224 I=ltNEVENT 
F2(I)=P1(I) 

P2( I)=P1( I) 

P12( I)=P11( I) 

F12( I )=P11 ( I) 

P = (I)=P4( I ) 

2224 P6( I )=P4( I ) 

GCTO 1111 

C 

C CALCULATE P20 AND P50 AND P120 VECTORS 

C 

2225 PF=1. 

N=NEVENT 

PFF=FLOAT (N+1) /FLOAT (N+2) 

P = l. 

J = C 

DC 2399 1=1 , NEVENT 

IF ( IT( I ) .EQ.l ) GCTO 2350 

J = 1 

GCTO 2399 

235C PP = PP=!'FL0AT{N-1) /FLGAT(N) 

P = P=)=FLOAT{N )/FLOAT(N + l) 

PFF =F FP* FLOAT! N)/FL0 AT (N+1) 

P2( I )=PP 
F5( I )=P 
P12(I)=PPP 
J=C 
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2399 N=N-1 

IF (J.EQ.O) GOTC2A25 
IF (NOI .GT.l) GOTO 2400 
T;flL=-TI I/ALCG( P5 ( II ) ) 

TTAU=-TII/ALOG( P12 (II ) ) 

GCTC241 0 

24 CC TAU = DTI /A LOG ( P3 ( I I ) / P5 ( 1 1 I) ) 

TTAU=DTI/ALCG(P12( II )/P12( I II ) ) 

241C IF (0T.GT.150=(^TAU) TAU = 0T/150 
P5(NEVENT) = F=i'EXP (-DT/TAU) 

I F(OT.GT. 15C-TTAU )TTAU = DT/i 50. 
P12(NEVENT)=FPP^EXP(-CT/TTAU ) 

C 

C CALCULATE P3() ANC P60 ANO PIBO VECTORS 

2425 PF=1. 

N=NEVENT 
P = l. 

PFP=FLOAT(N+l) /FLCAT(N+2) 

J = C 
TT = 0 
TTT = 0 

DC 2599 I=1,NEVENT 

IF ( IT( I ).EG.l ) GCTO 2550 

J = J + 1 

TT = TT+T( II-TTT 
GCTO 2599 
255C DK=N 

IF (TT.NE.O) ON=DN + TT/(T( I )-TTT ) 

PF = PP=!=( DN-1 )/0N 
P=F*DN/ (DN+1) 

FFP=PPP*DN/ (DN + 1 . ) 

P5( I )=PP 
P6( IJ=P 
P15( IJ=PPP 
J = C 
TT = 0 

TTT=T( I ) 

2599 N=N-1 

IF(J.EQ.O) GO TO 1111 
CTT=TN-TTT 

IF( CTT.GE.1E-70)GGTC 3005 
CN=.5=!'( J + 1 ) 

GCTC 2010 
30C5 DN=TT/DTT 

30 1C P6(NEVENT)=P*SQRT(CN/ (DN+1. ) ) 

IF(NOI.GT.l )G0 TO 1997 
TTAU=-TII/AL0G(P13(II ) ) 

GCTO 1998 

199 7 TTAL=DTI/ALCG (P13 ( II ) /P13 ( I II ) ) 

1998 IF(DT.GT.15C*TTAU) TTAU=DT/150. 

P13(NEVENT )=PPP=f=EXP(-DT/TTAU) 

C 

C SET UP A LOOP FOR ALL JACKKNIFE CALLS ( P J4 , P J5 , F J6 ) 
C 

llii CC 1000 I=1,NEVENT 
DC 1011 J=1,NEVENT 
PJ2( I , J ) = 0.C 
PJ4( I , J)=O.C 
Pj5(I, J)=O.C 
P J6( I , J J=O.C 
SLKI , J)=O.C 
SL2( I, J J = 0.G 
SL3( I , J)=O.C 
SL4( I , J ) = 0. C 
1011 CCKTINUE 

PJA2( I ) =0.0 
PJA4(I )=0.0 
PoA5(I)=0.0 
PJA6(I )=0.0 
SEARl ( I J=0.C 
SBAR2( I )=0.C 
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SE;iR3( I )=0.C 
S8AR4( I J =0 .C 
PZ2( I )=0.0 
VARJ2( I )=0.C 
VARJ4< n=o.c 
VARJ5( I }=0. C 
VARJ6( I )=0.C 
lOCC CONTINUE 

DC 1001 I=1,NEVENT 
C 

C MOVES DATA INTC TJO AND ITJO VECTCFS 
C 

K = 1 

JNEXT=0 

JEEFCR=0 

JAFTER=0 

10C2 IF(K.EQ.I )GC TO 7003 
TJ(K)=T (K ) 

I IJ(K)=IT(K) 

IF( IT(K } .EQ.OIGO TC 7001 
JNEXT=JBEFCP 
JEEFCR=K 
70C1 K=K+1 

GC TO 1002 
70C2 JAFTER=JBEFCR 
JEEFCR=JNEXT 
GC TO 1010 

70C2 IF(I.GT.JEVENT) GO TC 7002 
1002 IF(K.GT.JEVENT)GO TO 1004 
TJ(K)=T(K + 1 ) 

ITJ(K)=IT( K+1) 

IF(ITJ(K) .EC.OIGO TO 4002 
IF( JAFTER.EC.O ) JAFTER=K 
40C2 K=K+1 

GC TO 1003 

lOCA IF(JAFTER .E G .0 ) J A FT ER= JEVENT 
lOlC NDIJ=NDI-IT(I) 

C 

C CHECK IF ZERO DEATHS 
C 

I F (NCIJ .EQ.O )GOTO 1001 

N=NDI J 

P=l. 

J = C 
11 = 0 
I 11 = 0 

GC TC 1014 
C 

C CALCLLATE PJ4()VECT0RS 
C 

1014 DC lOlo IJ=1, JEVENT 
IF(ITJ( IJ).EQ.i)GG TO 1015 
J = 1 

GC TO 1016 

1015 P=FL0AT(N)/FL0AT(NCIJ+1) 

pzn J) =p 

I II=I I 
I I = IJ 
J = 0 
N=N-1 

1016 CONTINUE 
IF(J.EQ.O)GC TO 1019 
TII=TJ( II } 

CTI=TJ( III)-TI I 
I F(NDI J .GT.l JGO TC 1017 
TAU=-TI I/ALCG( PZ( ID) 

GC TO 1018 

1017 tal=oti/alog(p z( ii)/pzniin 

lOie DT=TJ( JEVENT)-TI I 

I F(DT.GT.150.-TAU )TAU=0T/150. 

PZ( JEVENT) =F*EXP(-0T/TAU) 

1019 K=1 
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x02C IF(K.EQ.I)GC TO 2021 

PJ^(K, I J = ALCG(( P2(K}+A)/(1-PZ(K ) + /i) ) 

K = K + 1 

GC TO 1C20 

2021 IF(K.6Q.NEVEMT)G0 TO 5021 
IF( IT(I ).EQ.O)GO TC 1025 
5021 TJA=TJ( JAFTER) 

IF(JBEFCR.NE.O)GO TO 1022 
PX = ALOG( P2( JAFTER) )*T(IJ/TJA 
PZZ=0.0 

IF (PX.LT.-150) GOTO 7025 
PZZ=EXP(PX) 

7025 PJ4(I,I1=ALCG( (PZZ+A)/(1-PZZ+A) ) 

GC TO 1025 
1022 T^B=TJ ( J8EFCR) 

DTJ=TJA-TJ8 

CTE=T(I)-Tje 

PX=DTB«ALOG (PZ( JAFTER )/PZ(J6EFCR))/0Tj 
PZZ=0.0 

IF (PX.LT.-150) GOTO 7026 
PZZ=PZ( JBEFCR)*EXP(PX) 

70 2 6 PJ4( I , I J = ALCG( ( PZZ + A)/( 1-PZZ+An 

1025 IF(K.GT.JEVENT)GOTC 2027 

P J4(K + 1 ,I )=ALOG ( ( FZ(K) + A)/ ( 1-PZ (K )+A) ) 
K = K+1 

GC TO 1025 

3027 IF(NDI.NE.NEVENT)GCTC 1026 
DC 3028 K=1,NEVENT 
CC 3028 L=i,NEVENT 
PJ5(K,L)=PJ4(K,L) 

302£ PJ6 (K,L ) = PJ4(K, L ) 

GOTO 1001 
C 

C CALCULATE PJ5() VECTORS 
C 

1026 N=JEVENT 
P = l. 

PP = 1. 

J = C 

DC 1028 IJ=1,JEVENT 

IF( ITJ( IJJ.EQ.DGC TO 1027 

J = 1 

GC TO lC28 

1027 P=P*FL0AT(N)/FL0AT(N+1) 

PF = PP=i<FLCAT (N-1 )/FLOAT(N) 

PZ( IJ)=P 

PZ2( IJ )=PP 
J=0 

102£ N=N-1 

IF(J.EQ.O)GC TO 1031 
I F(NOIJ .GT.l )G0 TC 1029 
TAL=-TI I/ALCG( PZ( ID) 

GC TO 1030 

1029 TAL=DTI /ALCG(PZ( II)/PZ(III)J 
10 2 C IF(DT.GT.150*TAU) TAU=DT/150. 

PZ( JEVENT )=F*EXP(-CT/TAU) 

1031 K=1 

1022 IF(K.EQ.I)GC TO 2033 

PJ5(K,I ) = ALCG((PZ(K)+A)/(1-PZ(K )+A)) 
PJ2(K,I )=ALCG( (PZ2(K)+A)/ (1-PZ2 (K)+A) ) 
K = K + 1 

GC TO 1032 

2032 IFIK.EQ.NEVENTJGO TO 5033 
IFdTU ). EQ.OGO TO 1136 
5033 IF( JBEFCR.NE.OIGO TO 1135 

PX = T( I ) =!'A LOG (PZ( J AFTER) )/TjA 
P2Z=0.0 

IF (PX.LT. -150) GCTO 7027 
PZZ=EXP (PX) 

70 2 7 PJ5 ( I , I ) = ALCG( ( PZ Z + A)/ (1-PZ Z+A) ) 

GC TO 1136 

1125 PX = DTB#ALOG(PZ( JAFTER )/PZ( JBEFCR) )/DTJ 



68 




I ) 




ooo 



PZZ=0.3 

IF (PX.LT.-150) GOTO 70Z8 
F 2Z = PZ( JEEFCR)=f'EXP(PX ) 

70 2 S PJ5( I ,I )=ALCG( ( PZZ+A) / ( 1-PZZ + An 
1136 IFIK.GT.JEVENT )G0 TO 1036 

PJ5(K+1 ,I )=AL0G( ( PZ(K J+A)/ I 1-PZ (K )4A) ) 

PJ2(K+1,I )=AL0G( (PZ2(K)+AJ/(1-PZ2(K)+A)) 

K = K + 1 

GC TO 1156 

CALCLLATE PJ6() VECTORS 

1036 N=JEVENT 
P = l. 

J=0 
TT = C 
TTT=0 

DC 1038 IJ=1,JEVENT 

IF( ITJ( I J) .EQ.l )GC TO 1037 

J = 1 

TT=TT+TJ( IJ )-TTT 
GC TO 1058 

1037 CN=N 

IF (TT.NE.O) CN=DN+TT/ ( TJ ( IJ )-TTT ) 

P=P*ON/ (ON+1.) 

PZ( IJ)=P 
J=C 
TT = C 

TTT=TJ( IJ J 

1038 K=N-1 

IF(J.EQ.O)GC TO 1041 
CT=TJ( JEVENT)-TTT 
IF(DT.GT.1E-70)GC TO 1059 
0N=.5«( J+1) 

GC TO 1040 
1C39 ON=TT/DT 

1C4C PZ( JEVENT)=F*SQRT(CN/(0N+1. ) ) 

1041 K=1 

1042 IFIK.EQ.I )GC TO 2043 
PJ6(K,IJ=ALCG((PZ(K)+A)/(1-PZ(K)+A)) 

K=K+1 

GC TO 1042 

2043 IF(K.EQ.NEVENT)GO TO 5043 
IF( IT(I ).EG.O)GO TC 1146 
5043 IF( JBEFCR.NE.OJGG TO 1045 

PX=T(I) «ALOG(PZ( JAFTER) )/TJA 
PZZ=0.0 

IF (PX.LT.-150) GOTO 7029 
FZZ=EXP (PX) 

7029 PJ6(I tl )=ALCG( ( PZZ+A)/(i-PZZ+A) ) 

GC TO 1146 

10 45 PX = DTB*ALOG(PZ( JAFTER )/PZ( JBEFOR) )/OTJ 
PZZ=0.0 

IF (PX.LT.-150) GOTO 7030 
PZZ=PZ( JBEFCR)*EXP(PX) 

7030 PJ6( I ,I )=ALCG( ( PZZ+A) /( 1-PZZ+A) J 
1146 I F(K.GT.jEVENT)GOTO 1001 

PJ6(K + 1 ,I )=ALOG( ( FZ(K )+A)/ (1-PZ (K )+A) ) 

K = K + 1 

GC TO 1146 
1001 CCNTINUE 

DC 8888 I=1,NEVENT 
DC 8811 J=1,NEVENT 
I F ( I .EQ .J ) GCTO 8811 

SLK I, J )=FLCAT( NE VENT)*ALOG( ( P4 ( I ) +A ) / ( 1 -F4 ( I ) +A) ) - FLC 
«AT(NEVENT-1 )*PJ4( I ,J) 

SL2( I , J )=FLCAT (NEVENT)*ALOG( (P5(I)+A)/(i-P5(I)+A))-FLC 
*AT (NEVENT-1 )*PJ5( I , J ) 

SL3( I , J)=FLCAT (NEVENT )-ALOG( (P6(I)+A)/(1-P6(I)+A))-FLC 
*AT(NEV5NT-1 )«PJ6( I ,J) 

SL4( I , J ) = FLCAT( NEVENT )>!'-ALOG ((P2(I)+A)/(i-P2(I) + A))-FLC 
*AT(NEVENT-i )^PJ2 ( I , J) 
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8811 CONTINUE 
88ce CCNTINLE 

DC 9999 I=1,NEVENT 
DC 9911 J=1,NEVENT 
S8/^R1(I )=SL1( I , J ) /NEVENT+SbARl (I ) 

SeAR2(I)=SL2(ItJ) /NEVENT+SBAR2( I ) 

S8AR3 ( I )=SL2( I , J ) /NEVENT+S 8AR3 ( I ) 

SEAR4( I )=SL4( I , J ) /NE VEMT+S B AR4 ( I ) 

9911 CCNTINUE 

IF(SBAR1( !) .LT.-ISO. ) SBARK I)=-1£C. 

IF( SBARK I) .GT. 174.) SBARK I ) = 174. 

I F(SBAR2( I) .LT.-180. ) SBAR 2 ( I ) =- 1 fiC . 

IF(SBAR2(I) .GT.174.)SeAR2(I )=174. 

IF(SBAR2( I) .LT.-ieC. ) SBAR5( I)=-18C. 

IF(SBAR3 ( I ) .GT .174.)SEAR3( I ) = 174. 

IF(S8AR4( I) .LT.-l 80. )SBAR4( I J=-180. 

IF(SBAR4( I) .GT.174.)SBAR4(I ) = 174. 

EL1=EXP(SBAR1 ( I ) ) 

EL2=EXP( SBAR2( I ) ) 

EL3 = EXP(SBAR3( I ) ) 

EL4=EXP( S6AR4( I ) ) 

PJA4(I)=EL1/(1.+EL1) 

PJA5(I)=EL2/(1.+EL2) 

PJA6(I)=EL3/(i.+EL3) 

PJA2 (I ) =EL4/( 1 . + EL4) 

999= CCNTINUE 

DC 8999 I=1,NEVENT 
DC 8911 J=1,NEVENT 
IF(SL4( I, J) .EQ.O.OIGOTO 5111 

VARJ2( I ) = ( SL4( I , J )-SBAR4( I ) ) ** 2/ FLC AT ( NE V ENT-1 ) +V A R J2 ( 
^ I ) 

3iil"lF(SLl(I , J) .EQ.O.OIGOTO 3112 

V/5RJ4(I ) = (SL1( I , J)-SBARKI ) ) #*2 / FLCAT ( NE V £N T- 1 ) +'/ AR J4 ( 
^I ) 

3112 IF(SL2( I, J) .EQ.O.OGOTO 3115 

VARJ5(I)=(SL2( I,J )-S8AR2(I) ) =i=*2 / FLC AT ( N£ V ENT-1 )+ VAR J5 ( 

n' I } 

31i3"lF(SL5( I ,J) .EQ.O.OGOTO 8911 

VARJ6(I ) = (SL3( I, J )-S8AR3(I) )**2/FLCAT(NEV5NT-l)+VAPJ6( 
) 

8911 CCNTINUE 

RPS1(I)=SQRT(VARJ2(I) /FLOAT ( NEVEN 1 ) ) 

RNS2 ( I ) =SQRT(VARJ4( I ) /FLOAT ( NE VENT ) ) 

RNS3( I ) =SQRT(VARJ5( I ) /FLOAT ( NEV ENT ) ) 

Ry.S4( I ) =SQRT (VARJ6( I ) / FLOAT ( NEV ENT) ) 

UP1(I)=SBAR4(I)+RNS1(I)*TN 

ROKI )=SBAR4(I)-RySl(I)*TN 

UF2(I)=SBARl(I)+RyS2(I)-TN 

RC2(I ) =SBARi (I)-RyS2 ( I )*TN 

UF3( I ) =SBAR2( I)+pyS3( I )«TN 

RC5(I ) = SBAR2( I)-RyS3( I)*TN 

UP4( I)=SBAR3(I)+RyS4(I)-TN 

RC4(I) = SBAR;{I)-RNS4( D^i'TN 

IF(UP1( I ) .GT.174. )UP1( I )=174. 

IF(UP1( I ) .LT.-IOO. )UP1(I)=-100. 

IF(R01( I ) .GT.174. )R01( I )=174. 

IF(ROK I) .LT.-IOO. )R01(I)=-100. 

IF (UP2( I ) .GT. 174. )UP2( I )=174. 

IF(UP2( I ) .LT.-IOO. )UP2( I )=-100. 

IF(R02( I) .GT.174. )R02(I)=174. 

IF(R02( I) .LT.-IOO. )R02( I )=-100. 
IF(UP3(I).GT.174.)UP3(I)=174. 

IF(UP3( I ) .LT.-IOO. )UP3(I)=- 100. 

IF(R03( I ) .GT.174. )R03( I )=174. 

IF(R03( I ) .LT.-IOO. )R03 ( I )=-100. 

IF(UP4( I).GT.174.)UP4(I )=174. 

IF (UP4( I ) .LT.-IOO. )UP4( I )=-100. 

IF(R04( I) .GT.174. )R04(I) = 174. 

IF(R04( I ) .LT.-IOO. )R0 4( 1)=- 100. 

UINTl ( I )=1./(1.+1 ./EXP(UP1( I ) ) ) 
RINTl(I)=l./(l.+l./EXP(R01(I))) 

UINT2(I ) = 1./(1.+1 ./EXP(UP2( D) ) 
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RINT2( I )=1./{1. + 1./EXP(RQ2(I) )) 

UINT3( I ) = 1./(1. + 1/EXP (UP3(I ) ) ) 

RIM3(Ii = l./(l.+l ./EXP (R03 ( I ) ) ) 

UINT4( I)=1./{1.+1./EXP(UP4(I))) 

RINT4( I ) = 1./(1 .+1 ./EXP(R04( I ) ) ) 

89^9 CCNTIN'JE 
C 

C PRINT OUTPUT 

C 

3850 IF (NWRITE.EQ.OJ GOTO 3500 

UPIT£(6,8)(IREPL,IT(IJ,T(I),Pl(I)tP2(I)tP3(I)»P4(I),P5 
n »P6( n , PJA2( I ) ,PJ/^4( I ) , PJA5 ( I ) ,PJA6( I ) ,P11 ( I ) t F12( I 
tP13(I )tI = ltNEVENT) 

WRITE ( 6tlOi 
C 

C RECLCE VECTORS 

C 

350C K=1 

CC 4000 I=1»NEVENT 

IF ( IT( I J .EC.l ) GCTC 3900 

IF ( I .NE.NEVENT) GOTO 4000 

39CC T(K)=T( I) 

Pi(K)=Pl ( I ) 

P2(K)=P2( I ) 

F3(K)=P3( I ) 

PA( K) =P4( I ) 

PS(K)=P5( I ) 

P6(K)=P6( I ) 

FJA2( K) =PJA2(I ) 

FoA4( K) =PJA4( I ) 

PwA5 (K) =PJA5( I ) 

P JA6( K) =PJA6 (I ) 

P11(K)=P11 ( I ) 

P12(K)=P12( I) 

P13(K)=P13( I ) 

LINTl ( K )=UINT1 ( I ) 

RINTl (K)=RINT1 ( I ) 

LINT2(K )=UINT2( I ) 

RINT2 (K )=RINT2 ( I ) 

UNT3(K) =UI NTS ( I ) 

RINT3(K )=RINT3( I ) 

UINT4(K )=UINT4( I ) 

RINT4(K)=RINT4( I ) 

K = K + 1 

40CC CONTINUE 
C 

C CALCULATE OIFFERENCESt MEAN, RMSt AND MEAN ABS 

C 

K = 1 
TTT=0 
CT=-T(1 ) 

PP1 = 1 
PF2 = 1 
FF3=1 
PFA = 1 
FF5=1 
PF6=1 
PF8 = 1 
PF9=1 
PF1C=1 

PF11= (NDI + l . )/ DI I 

PF12= (NEVENT+1. )/ (NEVENT+2. ) 

PF13=PP12 

FF14=1. 

PG4=P4( 1) 

PC5=P5( 1) 

FC6=P6( 1) 

FG8=PJA4 (1 ) 

P09=PJA5( 1 ) 

PG10=PJA6 (1 ) 

CCNF1 = 1 . 

CCNF2=1 . 
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CCNF3=1 . 

CCKF4=1 . 

CCNF5 = 1 . 

CCNF6=1 . 

CCNF7=1 . 

CCNF8=1 . 

7^=DT/AL0G( FC4) 

T5=DT/AL0G( PG5 ) 

T6=DT/AL0G( FC6 ) 

T£=DT/ALOG( FG8) 

T9=DT/AL0G(PQ9) 

T1C=DT/ALGG(PQ10 ) 

DC 4998 J=l,9 

TT=.1*J 

P=1-TT 

IF (NDIE.EQ.2) TT=-ALOG(P) 

IF (TT.GT.TN) GOTC 4999 
NN ( J)=NN( J)+l 

4i00 IF (TT.LE.T(K)) GCTC 4200 
PF1=P1 ( K) 

FF2=P2( K) 

FF2=P3( K) 

FF4=PQ4 

PP5=PQ5 

FF6=PQ6 

FF8=PQ3 

PF9=PQ9 

PF10=PQ10 

PFll=Pil(K) 

PF12=P12( K) 

PF13=;P12(K) 

PF14 = P JA2 (K ) 

CCNF1 = UINT1 (K) 

CCNF2=R INTI (K) 

CCNF3=UINT2 (K) 

CCNF4=R INT2(K) 

CCNF5=UINT3 (K) 

CCNF6=RINT3 (K) 

CCNF7=UINT4(K) 

CCNF8=R INT4 (K) 

TTT=T(K) 

K = K + 1 
PC4=P4( K) 

PC5=P5( K) 

PC6=P6( K) 

PG8=PJA4(K) 

PC9=PJA5(K) 

PC10=PJA6(K) 

D7=T(KI -TTT 
TT1=AL0G( PP£/PQ8) 

TT2=AL0G( PP9/PQ9 ) 

T73=ALOG( PFIO/PQI C ) 

IF(PQ8.EQ.O.C.OR. TT1.EQ.0.0J7T1=ALCG( ( PP8+A) /(1-PG8+A)) 
IF(PC9. Ey.O.O.OR.TT2.EQ.O-0)TT2=ALCG( ( PP9+ A ) / ( 1-PQ9 + A )) 

IF (PgiO.EQ.C.O.OR.TT3.EQ.O.O) TT3 = ALaG( (PP10 + A)/(l-PglC+A)) 
T4=DT/ALCG( PP4/PQ4) 

T5=DT/ALCG( FF5/PQ5 ) 

T£=DT/ALCG( FP6/PQ6) 

T£=DT/TT1 
79=07/772 
T1C=DT/TT3 
GCTC 4100 

42CC D(1)=FP1-P 
C ( 2) = PP2-P 
C (3) = FP3-P 
017=777-77 

0 (4 ) = PP4*EXF( DTT/T4) -P 
C(5)=PP5=!=EXF(DTT/T5)-P 
D( 6 )=PP6*EXF( 077/76) -P 
D(7) = EXP(TT=^T7)-P 
C(£J=PP8*EXF(DT7/T3)-P 
D( 9) = PP9'-!'EXF{DTT/79)-P 
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C (10)=PP10*EXP( DTT/TIO )-P 

D( 11 )=PP11-F 

D (12)=PP1Z-F 

C(13)=PP13-F 

0 ( 14)=PP14-P 

CF (1 )=CCNF1 

CF (2) =CCNF2 

CF (3J=CCNF3 

CF(4)=CCNF4 

CF (5)=CCNF5 

CF (6)=CCNF6 

CF (7)=CCNF7 

CF( 8J=CCNF8 

DC 491 1 1=1,8 

CC=CF( I ) 

C ( I ,J)=C( I , JJ+CC 
IF(NWRITE.EG.O)GQTC 4911 
UPITE(6 ,10) 

'wRITE(6,12)I,P,(C(I,L),L=l,9) 

4911 CCMINUE 

DC 4997 1=1, NP 

DC = D( I ) 

U(I,JJ=U( I, J)+DD 
V>(I,J)=^^( I, J)+ABS(CD) 

S(I,J) = S( I, J)+OD=!‘DD 
IF (NWRITE. EQ.O ) GCTO 4997 
VvRITE ( 6,10 ) 

WRITE (6,12) I,P ,(U(I,L) ,L=1,9) 

WRITE (6,12) I,P , (W( I ,L) ,L = 1,9 ) 

WRITE (6,12) I,P ,(S(I,L) ,L=1,9) 

4997 CONTINUE 

4998 CCNTINUE 
C 

4999 CCNTINUE 
C 

C PRINT SUMMARY STATISTICS 

C 

DC 5100 J=1 ,9 
51CC 0(J)=.1*J 

WRITE (6,7) (D( J) ,J = 1,9) 

WRITE (6,10 
DC 5300 J=l,9 
IF (NN(J).GT.O) GCTO 5150 
DC 5125 1=1 ,NP 
S( I,J) = 1E9 
U(I,J)=1E9 
C( I,J)=1E9 
5125 W( I,J)=1E9 
GCTC 5300 
515C XJ=1./NN(J) 

52C0 DC 5250 1=1, NP 

S(I,J)=SGRT(S( I, J)^XJ) 

U(I,J)=U(I, J)*XJ 
C(I,J) = C( I, J)*XJ 
525C W(I,J)=W( I,J)TXJ 
52CC CCNTINUE 

WRITE (6,9) (I, (U(I, J), J=l,9) ,I=1,NP) 

WRITE (6,10) 

WRITE (6,14)( I, (W(I, J),J=1,9) ,I=1,NP) 
WRITE (6,10) 

WRITE ( 6,16)(I, (S(I,J), J=l,9) ,I=1,NP) 
WRITE (6,10) 

WRITE(6,20) (I, (C(I ,J) ,J=1 ,9 ) ,1 = 1,8 ) 

20 FCRMAT( IX, • C ,11, 9F8.3,3X, • CCNF* ) 
wRITE(6 ,10 ) 

WRITE (6,15) (NN( J) , J=1 ,9) 

IF (NENC.NE.O) GOTC 108 
WRITE ( 6,11 ) 

GCTO 25 
ENC 
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